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Abstract. In the past years we have experienced an increasing interest in 
understanding of the physical properties of collisionless plasmas, mostly because of 
the large number of astrophysical environments, e.g. the intracluster medium (ICM), 
containing magnetic fields which are strong enough to be coupled with the ionized gas 
and characterized by densities sufficiently low to prevent the pressure isotropization 
with respect to the magnetic line direction. Under these conditions a new class of 
kinetic instabilities arises, such as firehose and mirror ones, which were extensively 
studied in the literature. Their role in the turbulence evolution and cascade process in 
the presence of pressure anisotropy, however, is still unclear. In this work we present the 
first statistical analysis of turbulence in collisionless plasmas using three dimensional 
double isothermal magnetohydrodynamical with the Chew-Goldberger-Low closure 
(CGL-MHD) numerical simulations. We study models with different initial conditions 
to account for the firehose and mirror instabilities and to obtain different turbulent 
regimes. We found that the CGL-MHD subsonic and supersonic turbulence show 
small differences comparing to the MHD models in most of the cases. However, in the 
regimes of strong kinetic instabilities the statistics, i.e., the probability distribution 
functions (PDF) of density and velocity are very different. In subsonic models the 
instabilities cause an increase in the dispersion of density, while the dispersion of 
velocity is increased by a large factor in some cases. Moreover, the spectra of density 
and velocity show increased power at small scales explained by the high growth rate of 
the instabilities. Finally, we calculated the structure functions of velocity and density 
fluctuations in the local reference frame defined by the direction of magnetic lines. The 
results indicate that in some cases the instabilities significantly increase the anisotropy 
of fluctuations. These results, even though preliminary and restricted to very specific 
conditions, show that the physical properties of turbulence in collisionless plasmas, as 
those found in the ICM, may be very different from what has been largely believed. 
Implications can range from interchange of energies to cosmic rays acceleration. 
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1. Introduction 

Magnetized and low density (weakly coUisional) plasmas are known to present 
anisotropic pressures with respect to the magnetic field orientation, which can survive 
considerably long times compared to the dynamical timescales of certain systems. 
In astrophysical environments such pressure anisotropy may be generated by several 
different processes, such as kinetic pressure of cosmic rays, supernovae explosions, stellar 
winds, or anisotropic turbulent motions (see Quest and Shapiro 1996). 

Under certain conditions gyrotropic plasmas give rise to new wave modes and 
instabilities, which cannot be studied by standard isotropic magnetohydrodynamic 
(MHD) model (Hascgawa 1969, Wang and Hau 2003, Passot and Sulem 2006). 
For instance, Hau and Wang (2007) showed that gyrotropic magnetohydrodynamic 
equations closed by the Chew-Goldberger-Low laws (CGL-MHD) lead to a positive 
density n versus magnetic field strength B correlations for the slow magnetosonic mode 
under certain conditions in contradiction to the standard MHD model. It was commonly 
believed that the detected absence of slow modes in several astrophysical sites was 
related to the strong damping of these waves. However, in the collisionless plasmas it 
could be explained by wrong identification of the positive n — B correlations for fast 
magnetosonic modes. 

The pressure anisotropics give rise to plasma instabilities depending on the 
anisotropy ratio, e.g., p\\ > p± and p± > for the firehose and mirror instabilities, 
respectively. They are responsible for the growth of the magnetic energy and acceleration 
of particles. The predictions of the CGL-MHD model, including new plasma instabihties, 
are also important in weakly magnetized environments, since even a weak magnetic 
field is enough to change the motion of the charged particles and therefore increase the 
pressure anisotropy. 

As an example of the possible applications Sharma et al. (2003) and Sharma et al. 
(2006) showed the importance of the collisionless plasma approach in protostellar disks. 
They studied the role of the kinetic instabilities in the magneto-rotational instability 
(MRl) showing that the transport of angular momentum in the disk may be efficiently 
increased under certain circumstances. 

The intracluster medium (ICM) is possibly the most suitable environment for 
studies of gyrotropic plasma effects (Schekochihin et al. 2005). Considering typical 
parameters of n ~ 10~ cm" -3, T ~ 10^ K and S ~ 1 ^G (Ensshn and Vogt 2006), it 
is possible to show that the cyclotron frequency {VL) is much larger than the collision 
frequency {vu). Under such conditions, the plasma fluctuations with wave numbers 
k >1 kpc will be subject to different processes related to the pressure anisotropy. 
The turbulent cascade, for example, may be modified by the new wave modes and 
instabilities resulting in a different picture of the energy budget in these environments. 
It may be particularly important in understanding of the cooling flow and the cosmic 
rays acceleration processes. 

Schekochihin et al. (2008) showed that in the ICM the anisotropy-induced flrehose 
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and mirror instabilities may grow non-linearly up to the saturation at SB/B ~ 1. 
This growth is very fast. They estimated the increase of power at small scales. As a 
consequence, the excess of small scale fluctuations of B and the energy transport at 
these environments may drastically change the picture. In certain cases the small scales 
may also be considered as collisionless, what may be important, for instance, in the 
development of turbulent cascade. For wavelengths smaller than the Larmor radius the 
kinetic treatment of plasma is necessary, as shown by Howes et al. (2006). However, 
in this work we focus on the large scale so the CGL-MHD approximation may be used 
instead. 

The aim of this work is to provide an extensive statistical analysis of the MHD 
turbulence in collisionless plasmas and to study the role of the different instabilities 
in the evolution of the system. To accomplish that we performed the first 3- 
dimensional simulations focusing on the evolution of turbulence in the presence of 
pressure anisotropy. The description of the model as well as the presentation of the 
governing equations is given in Section |2] with additional discussion of the instabilities 
and the double-isothermal approximation. In Section |3] we describe the numerical 
simulations. The results and the extensive statistical analysis including the derivation of 
probability distribution functions, spectra and structure functions of the fluctuations of 
density and velocity are presented in Section HI We discuss the most important results 
in Section [5] followed by Section |6l where we draw main conclusions. 



2. The Double-Isothermal CGL-MHD Approximation 



2.1. Governing Equations 

In the fluid approximation, a gyrotropic plasma can be described by the Chew et 
al. (1956) magnetohydrodynamic (CGL-MHD henceforth) equations expressed in the 
conservative form as follows: 
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where p and v are the plasma density and velocity, respectively, B is the magnetic field, 
P = p±i + (p|| — p±)bb is the pressure tensor, b = B/|B| is the unit vector along the 
magnetic field, p\\ and p± are the pressure components parallel and perpendicular to b, 
respectively, and / represents the forcing term. 

The above set of equations is completed by the description of the parallel and 
perpendicular pressure components. To avoid the complexities related to explicitly 
calculated processes that may be important for the description of the energies, e.g., 
the anisotropic heat conduction and the emission and absorption of the radiation, we 
can make use of the double-polytropic equations instead, as suggested by Chew et al. 
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(1956). Tests of this approximation using the solar magnetosheath data confirmed its 
validity under the form (see Hau and Sonnerup 1993): 

d ( p± \ 



dt \pB^^-^ J 
d fpiiB^w-^' 



0, (2a) 
0, (2b) 



dt \ p^ii 

where 7^ and 7j| are the polytropic exponents for the perpendicular and parallel 
pressures, respectively. These equations, according to Hau (2002), can be expressed 
in the conservative form, as follows: 

^ + V • (S^v) = 0, (3a) 
aq 

^ + v-(5||^) =0, m 

where Sj_ = p±B^~'^^, S\\ = p\\ {B/py^^~ , and -B = |-B| is the strength of magnetic field. 
In this form, the double-polytropic equations can be solved numerically similarly as the 
continuity equation. 

2.2. Double-Isothermal Closure 

Under the double-isothermal closure we have 7|| = 7^ = 1 and the equation of state 
is described by two relations, p± = a\p and p\\ = ayp, where a± and ay are constants 
and represent speeds of sound along the perpendicular and parallel directions to the 
magnetic field, respectively. In such situation the conservative form of the momentum 
equation for the double-isothermal CGL-MHD model can be rewritten as follows, 
dpv 



dt 



+ V 



alp + ^]i-{l-a)BB 



I (4) 



where a = ^ \J5\^ — /S^j = (C ~ 1)) C = is the pressure anisotropy ratio, and 

^11 and are the plasma betas in the parallel and perpendicular directions to the local 
field, respectively. 

The main consequence of the double-isothermal closure is that the pressure 
anisotropy is kept constant, i.e. ^ = const, independently of the evolution of the kinetic 
instabilities that may arise. Therefore, the stability condition is fulfilled by the local 
decrease of density and the increase of magnetic pressure and not due to the local 
changes of the pressure tensor. 

Furthermore, this work is focused on studying the differences of the turbulent 
regimes in the collisionless plasmas and the standard MHD turbulence using as the 
reference the studies done by Kowal et al. (2007), Kowal and Lazarian (2010), where 
the authors presented an extensive statistical analysis of density, velocity and magnetic 
field distributions in the isothermal MHD simulations. For this reason we chose the 
double-isothermal CGL-MHD model as the natural extension of the isothermal MHD 
in order to understand the importance of pressure anisotropy and its consequences in 
the turbulent plasmas. 
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Figure 1. Stability conditions for models with the sound speed ratio a^^/a± = 0.5 
(left panel) and a||/a^ = 2.0 (right panel). On the horizontal axis we show the plasma 



beta parameter related to the parallel pressure and defined as = Pmag/p\\ = ^c\/ap 
where Pmag = 5 I^P is the magnetic pressure. On the vertical axis we show the squared 
sinus of the angle between the directions of the perturbation and mean field. 



2.3. Firehose and Mirror Instabilities 

While in the MHD model the dispersion relations do not exhibit any instabilities, in the 
CGL-MHD equations the term of the pressure anisotropy introduces significant changes 
to the dynamical system. The detailed linearization and wave analysis of the double- 
isothermal CGL-MHD equations is provided in the Appendix. Here we describe the final 
dispersion formulas resulting from the presence of pressure anisotropy. For example, in 
the case of the incompressible Alfven mode, the dispersion relation can be written as 

{ioye)^ = cl{l~a)cos^e, (5) 

where a = ^ {l3\\ - I3±) and 9 is the angle between the mean field and the wave vector of 
the perturbation. Equation becomes negative when — > 2, what results in the 
occurrence of the firehose instability and the growth of the magnetic field fluctuations. 
As a consequence, the field hues bend and the magnetic pressure increases reducing the 
parameter a to its saturation value a ^ 2. This is known as the kinetic Alfven firehose 
instability. Two other instabilities related to the slow mode are called the compressible 
firehose and mirror instabilities and occur when p\\ > p± or p± > p\\, respectively, and 
the dispersion relation for magnetosonic waves becomes negative (see Appendix). 

In Figure [1] we show the stability regimes for two cases of the pressure anisotropy 
studied in this paper, a\\/aj_ = 0.5 in the left panel and a\\/a± = 2.0 in the right one 
(corresponding to p\\/p± = 0.25 and p\\/p± = 4.0, respectively), plotted as functions of 
the plasma beta (3\\ related to the parallel pressure and the squared sinus of the angle 
between the directions of the perturbation propagation and magnetic field. When the 
perpendicular pressure dominates (left panel) only the mirror instability corresponding 
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to the slow mode can occur. When perturbations propagate in directions almost 
perpendicular to the magnetic field the range of unstable grows. On the contrary, 
when the direction of propagation is close to the direction of magnetic field the instability 
does not occur. In the case of dominating parallel pressure (right panel) the plasma 
becomes unstable due to two firehose instabilities related to the incompressible Alfven 
and compressible slow modes, but the Alfven mode firehose instability extends over 
larger region of the parameter space up to (3\\ = 0.375 independently of the angle between 
the perturbation propagation and magnetic field directions. 

Both instabilities have the growth rate 7 = lm((X') larger at smaller scales, i.e., 
increasing with the wave number k. This dependence introduces a stiffness in the 
numerical integration since the micro instabilities, which may be due to the numerical 
imprecision, tend to grow very fast and destroy the configuration of the studied problem. 
As a consequence, the pressure anisotropy tends to disappear quickly and the problem 
of interest cannot be studied. Sharma et al. (2006) studying the magneto-rotational 
instability (MRI) in protostellar disks pointed out this important problem. They 
bypassed it by implementing a quasi-stability condition for the computational cells where 
the pressure anisotropy was larger than the threshold for instability. Physically, it can 
be understood as an almost "instantaneous" evolution of the system to the quasi-stable 
condition. Under such condition all anisotropic effects are kept at the large scales and 
the problem of interest could be analyzed. The downside of this method is that it 
artificially removes the free energy of pressure anisotropy without, as a counterpart, 
increasing the kinetic and/or magnetic energies. That is because, in the physical sense, 
the instability increases the magnetic energy and/or accelerate the gas. 

In the case of turbulence this kind of artificial removal of the energy at small 
scales may not be the best approach since the small scale structures in turbulence are 
generated by the cascade process operating in the turbulent models and therefore by 
artificially infiuencing the dissipation of energy we can distort its physical meaning and 
obtain incorrect conclusions. Nevertheless, under the double-isothermal approximation 
the pressure anisotropy is kept even after the growth of the small scale perturbations 
saturates, so the stability condition at large scales may lay in the unstable region even 
after the saturation at small scales. Therefore, under certain conditions, the double- 
isothermal approximation is valid and might be an interesting area of studies of the 
evolution of turbulence in a collisionless plasma. 

3. Numerical Simulations 

The simulations of turbulence in collisionless plasma were performed solving the set of 
double-isothermal CGL-MHD equations in a conservative form given by equations ( fTa|) - 
( ITcj) with an addition term / in the motion equation representing the turbulence driving. 
The numerical integration of the system evolution governed by the CGL-MHD equations 
were performed using the second order shock-capturing Godunov-scheme code (Kowal 
et al. 2007, 2009, Kowal and Lazarian 2010). We incorporated the field interpolated 
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Table 1. Description of the performed simulations of the double-isothermal CGL- 
MHD turbulence with the resolution 512'^. Initial density for all models was set to 
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constrained transport (CT) scheme (see, e.g., Toth, 2000) into the integration of the 
induction equation to maintain the V--B = constraint numerically, the general Harten- 
Lax-van Leer Riemann solver (Einfeldt et al. 1991) to obtain the numerical fluxes. The 
time integration was done with the second order Runge-Kutta method (see e.g. Press 
et al. 1992). On the right-hand side, the source term / represents a random solenoidal 
large-scale driving force. The rms velocity Sv is maintained to be approximately unity, 
so that V can be viewed as the velocity measured in units of the rms velocity of the 
system and B/ ^/Airp as the Alfven velocity in the same units. The time t is in units 
of the large eddy turnover time (~ L/5v) and the length in units of L, the scale of 
the energy injection. The magnetic field consists of the uniform background field and 
a fluctuating field: B = Bq + b. Initially, b = 0.0. We use units in which the Alfven 
speed ca = Bq/ y/A-Kp =1.0 and p = 1.0 initially. 

For our calculations, similar to our earlier studies (Kowal et al. 2007, Kowal 
and Lazarian 2010), the sound speeds and the strength of the external field Bq are 
the controlling parameters defining the sonic Mach number M.s = (Sv/a) and the 
Alfvenic Mach number A4a = {^v/ca), respectively. The angle brackets () signify 
the averaging over the volume. A^s < 1 and A^^, > 1 define subsonic and supersonic 
regimes, respectively, and J^a < 1 and Ma > 1 define another two regimes, sub- 
Alfvenic and super- Alfvenic, respectively. Since these two parameters are independent 
we can analyze, e.g., supersonic sub- Alfvenic turbulence, which signifies that A^s > 1 
and MIa < 1. In the case of the CGL-MHD simulations the regimes are defined by two 
sonic Mach numbers corresponding to the parallel and perpendicular sounds speeds. 

We drove turbulence solenoidally at wave scale k equal to about 2.5 (2.5 times 
smaller than the size of the box). This scale defines the injection scale in our models. 
We did not set the viscosity and diffusion explicitly in our models. The scale at which 
the dissipation starts to act is defined by the numerical diffusivity of the scheme. 

We performed six three-dimensional CGL-MHD simulations using the resolution 
512^ for different initial conditions, as shown in Tabled! We simulated the clouds up to 
^max ~ 6, i.e. 6 times longer than the dynamical timescale, to ensure a full development of 
the turbulent cascade. The computational time required for each CGL-MHD simulation 
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with intermediate resolution was equivalent to a MHD 512^, as performed by Kowal et 
al. (2007) and Kowal and Lazarian (2010). 

Models 1 and 2 are examples of weak turbulence and belong to subsonic and 
subAlfvenic regimes. The difference between these models is the pressure anisotropy 
accounting for the mirror and firehose instabilities, the comparison of both models gives 
us an insight into the different evolution of turbulence in these two cases. Similarly 
we calculated Models 3 and 4 for strong turbulence resulting in the supersonic and 
superAlfvenic regimes. Again, we tested these regimes for different instabilities. Model 
5 belongs to the subsonic and superAlfvenic regime representing the physical conditions 
of the ICM, and can be used as reference for further studies on this subject. Finally, 
Model 6 belongs to the supersonic and subAlfvenic turbulent regime and is the only 
simulation initiated with magnetic pressure larger than the thermal one. In this model 
all cells are initially stable but as the turbulence develops unstable conditions can occur. 
In the following sections we present the results obtained for each model, as well as a 
direct comparison with the standard MHD turbulent simulations with similar initial 
conditions. 

4. Results 

4.I. Distribution of Density and Velocity 

The results obtained for the distribution of density and velocity show strong differences 
between the CGL-MHD and standard MHD models. However, the kinetic instabilities 
play a role in the evolution of turbulence under specific conditions. In Figure |2] we 
present the column density obtained for each CGL-MHD model, as as well as the MHD 
models for comparison. Each case is presented as labeled in Table 1, with the MHD 
case shown in the left, and the CGL-MHD in the right. 

From the column density maps it is possible to distinguish models 2 and 5. Both 
models are subsonic, where initially we set p\\/p± = 2, which resulted in strong firehose 
instabilities. Here, the firehose instability is responsible for a deformation of the 
magnetic field lines. The curved magnetic lines tend to slow down and trap the flowing 
gas. Since the growth rate is larger at small scales we expect this effect to create more 
granulated maps. This is not seen in model 4 because pSv"^ > p\\ SB^ and, therefore, 
the turbulence is able to destroy the configuration of the growing instability. The same 
occurs for model 3. Even though it is not visible in the column density maps, the kinetic 
instabilities are responsible for changes in the statistics of the turbulence, as addressed 
below. For models 1 and 6, in which p\\/p± = 0.5, the column density maps show mild 
differences. In these cases we have a mirror instability operating, which is responsible 
for changes in the velocity distribution. For model 1 with a weak turbulence, the 
mirror instability is responsible for the acceleration of the gas resulting in the increase 
of the effective sonic Mach number. In model 6, where initially all cells were stable, 
the evolution of turbulence causes the instability in most of the computational domain. 
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The instability is responsible here for slowing the gas and reducing the effective sonic 
Mach number. 

In Figure [3] we present the probability distribution functions (PDFs) of density 
obtained from each model. In general, the distribution of gas exhibits an increasing 
contrast with increasing sonic Mach number. This result is in agreement with the 
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Figure 3. PDFs of density for tiie studied models. The left column shows models 1 
and 2, the middle one shows models 3 and 4, and the right column shows models 4 and 
5, according to Table [T] For each case both, the CGL-MHD (solid lines) and MHD 
(dashed lines) models are shown for comparison. 



Histogrom of Velocity 




Histogrom of Velocity 




Histogrom of Velocity 




Figure 4. PDFs of velocity for the studied models. The left column shows models 
1 and 2, the middle one shows models 3 and 4, and the right column shows models 
4 and 5, according to Table [TJ For each case both, the CGL-MHD (solid lines) and 
MHD (dashed lines) models are shown for comparison. 



MHD models. However, the kinetic instabilities cause even larger density contrast in 
the models with weak turbulence, compared to the MHD models. As stated above, 
in models 1 and 2 the instabilities could freely grow without being suppressed by the 
turbulent motions of the gas. The PDFs of density for the remaining models show no 
difference when compared to the MHD models. Surprisingly, in model 5 the firehose 
instability is responsible for the granulated map of column density, but the PDF of 



Turbulence in collisionless plasmas 



11 



density remains very similar to the MHD case. 

The PDFs of velocities are shown in Figure HI Again, the supersonic and super- 
Alfvenic models (3 and 4) show small differences compared to the standard MHD 
case. On the other hand, the subsonic and sub-Alfvenic models (1 and 2) present 
distorted velocity PDFs, with an increase in the high velocity tail. Both instabilities are 
responsible for more effective acceleration of the plasma, broadening the distribution of 
velocities. In the case of the firehose instability, the weakly magnetized model presents 
more evident changes, as noted in model 5. Here, since the magnetic field is weak the 
perpendicular flows are able to destroy the magnetic field configuration. Therefore, the 
magnetic breaking is not important, what results in flows with higher velocities. The 
same process is responsible for the changes in the PDF of model 6. Interestingly, for 
model 6 we obtained narrower PDF with lower velocities when compared to the MHD 
case. As the turbulence develops and unstable cells arise the local magnetic field grows. 
The strong magnetic breaking takes place resulting in a low velocity distribution. 

4-2. Power Spectra of Density and Velocity 

In order to characterize and study the changes in the energy cascade, as well as the 
correlations between different scales in CGL-MHD models, we analyze the power spectra. 
In Figures [5] and [6] we present the spectra of fluctuations of density and velocity for all 
models. Within each plot, we also present the anisotropy of the spectra for parallel and 
perpendicular directions with respect to the global magnetic field. 

At the center, we show the plots for Models 3 and 4. As also shown previously, 
there is no substantial difference between the CGL-MHD and the MHD simulations. 
The density spectra present similar slope {a ~ —0.5) at the inertial range, while a slope 
of ~ —2 is obtained for the velocity spectra. Regarding the anisotropy (/y oc /j^), we 
found the typical relations ( ^ 1 at small scales, but the Goldreich and Sidhar (1995) 
slope ( ~ 2/3 at larger scales, for both MHD and CGL-MHD simulations. A similar 
result in spectra and anisotropy of density was obtained for Model 6, with exception of 
a shght increase in the spectrum at small scales. 

Clearly, differences appear between the spectra of the MHD and CGL-MHD 
simulations in Models 1, 2 and 5. In these cases, it is noticeable the power excess 
at large values of k due to the fast growth of the instabilities in these scales. The 
slopes of velocity and density spectra for the MHD models range between -1.7 to - 
2.0 at the inertial range, while CGL-MHD simulations show positive slopes (~ +1 for 
Models 1 and 5) in velocity spectra, and fiat density spectra (a ~ 0). Regarding the 
spectral anisotropy, as an interesting result, we see that the firehose instability reduces 
the anisotropies regarding the global magnetic field lines. A more detailed study of the 
anisotropy of density and velocity perturbations regarding the local magnetic field is 
given below. 
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Figure 5. Spectra of density for tlie studied models, models 1 to 3 in the upper row 
and models 4 to 6 in the lower row, according to Table [T] For each case both, the 
MHD (solid lines) and CGL-MHD (dashed lines) models are shown for comparison. 



4-3. Structure functions 

In the previous section we showed, from power spectra of density and velocity 
fluctuation, that the instabilities in CGL-MHD models efficiently transport power from 
large to small scales, where their concentration is increased. Interesting modifications 
from MHD turbulence was also shown from the spectral anisotropy regarding the global 
magnetic field. However, since the instabilities are dominant at small scales, mapping 
the structure functions regarding the local magnetic field seems to be a better approach 
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Power Spectrum of Velocity: Be,,= 1.0, a||=1.0, aj_ = 2.0, model 1 
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Power Spectrum of Velocity: Bg,,= 0.1, O||=0.1, a^=0.05, model 4 
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Figure 6. Spectra of velocity for tlie studied models, models 1 to 3 in the upper row 
and models 4 to 6 in the lower row, according to Table [1] For each case both, the 
MHD (solid lines) and CGL-MHD (dashed lines) models are shown for comparison. 



if we want to determine the role of the instabihties on the isotropization of fluctuations. 
The second order structure function of a given parameter / is defined as: 

SF(0 = (|/(r + l)-/(r)|^), (6) 

where r is the referenced position and I is the distance calculated along the magnetic 
field line. The SF is calculated by randomly choosing a large number of referenced 
positions for each studied correlation length /. Here, to account for the importance of 
magnetic field fluctuation at small scales, we calculate the correlation length / along the 
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Figure 7. Structure functions of the density in the local reference frame for the studied 
models. In the left column we show models 1 and 2, inte middle column models 2 and 
3, and in the right column, models 5 and 6, according to Table [1] For each case both, 
the CGL-MHD (solid lines) and MHD (dashed lines) models are shown for comparison. 



The obtained SFs for density and velocity fluctuations are shown in Figures [7] and 
m respectively. In all plots, the solid lines represent the results obtained from the MHD 
simulations, while the dashed lines were obtained from the CGL-MHD simulations. The 
numbers refer to each model, as described in Table [TJ 

Similar to the previous calculations. Model 3 shows no difference between the two 
theoretical approaches, independent on the scale. The same result is obtained for Model 
6, which also showed similar spectra for the CGL-MHD and MHD models. Surprisingly, 
the structure functions of Model 4, which showed no differences in PDF and spectra, 
present more isotropic maps in the CGL-MHD model. The same behavior is obtained 
for Models 2 and 5. These are the models presenting p\\/p± = 2, i.e. the firehose 
instability. Here, the firehose instability is responsible for changes in the magnetic field 
topology, tangling the field lines, resulting in an increase of the perpendicular pressure 
in the local reference frame. Obviously, this effect is over estimated in these calculations 
because of the double-isothermal approximation. Otherwise, the magnetic field topology 
would not change drastically, but the interchange of energy would cause a reduction of 
parallel pressure anyway. In this sense, the result would be the same, i.e. the firehose 
instability is responsible for a isotropization of the fluctuations with respect to the 
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Strjcture Function: Velocity Structure Function; Velocity Structure Function: Velocity 

model 1: 3^,^ = 1-0, a,, = 1.0, Oi = 2.0 model 3; B„t = 0.1, 0, = 0.1, Oi = 0.2 model 5: B„j = 0.1, a„ = 1.0, o^ = 0.5 




Figure 8. Structure functions of velocity in the local reference frame for the studied 
models. In the left column we show models 1 and 2, inte middle column models 2 and 
3, and in the right column, models 5 and 6, according to Table [1] For each case both, 
the CGL-MHD (solid lines) and MHD (dashed lines) models are shown for comparison. 



magnetic field lines. Model 1, on the other hand, presented a larger anisotropy for the 
CGL-MHD model. Here, as p\\/p± = 0.5, the free energy is mostly converted to kinetic 
pressure. The mirror instability is responsible for an increase in the acceleration of the 
plasma along the field lines (as already seen in the PDF of velocity), increasing the 
anisotropy of the fluctuations. The result is more elongated structures, mostly at small 
and intermediate scales, as also noticed in the power spectra. 

5. Discussion 

5.1. Time Evolution of Unstable Systems 

In order to understand the evolution of the system during the growth and saturation 
of the instabilities, we calculated the instability condition for the Alfven and the 
compressible modes for each cell of the computational domain. The dispersion relation 
of these modes reveal the cells where the instability condition is fulfilled. 

As a general behavior, we found that the turbulence increases the range of unstable 
cells, even for the case with initially stable cells (Model 6). This process is illustrated in 
Figure [9] where we show the time evolution of the dispersion relation of Model 4 (similar 



Turbulence in collisionless plasmas 



16 



Firehose — Instability Condition Firehose — Instability Condition 




-20 -15 -10 -5 5 -0.04 -0.02 0.00 0.02 0.04 



Figure 9. Probability distribution function of stability condition for incompressible 
Alfvenic (left) and the slow magnetosonic (right) waves. 

distributions were obtained for all models). Here, we plot the histograms of the stability 
condition for different times. The Alfvcn and magnetosonic modes are independently 
shown in the left and right plots, respectively. Unstable cells populate the negative 
range of the dispersion relation. The time is shown in units of the dynamical timescale 
Td ~ 5V/L. 

As the turbulence is injected in the simulation domain, fluctuations of density and 
magnetic field change the local characteristic speeds and, as a consequence, the stability 
conditions. This process is fast compared to the dynamical timescale of the system 
{t < O.bTd). As a result, the dispersion relation spread over a larger range including 
the negative values of {u/kY. Then, the instabilities start to grow and saturate, at 
different timescales for different scales. The saturation of the instability brings the cells 
towards positive values of (w/Zc)^. At t ~ 2.5, most of cells have already reached the 
saturation condition for both modes. We believe that the few still unstable conditions 
are related to the large scale fluctuations, which evolve slow and has higher saturation 
values. It also seems that the firehose instability of Alfven modes saturates faster than 
the magnetosonic mode, though the differences between the modes could be also be 
related to the driving mechanism. We plan to address this possibility in a future work. 

In the present calculations we did not perform any wave mode decomposition and, 
therefore we were unable to analyze the time evolution of the individual modes. Also, 
for this reason we are unable to characterize the time evolution for each wavelength 
(scale). Needless to say that this is an interesting subject for future studies on CGL- 
MHD turbulence. 

6. Conclusions 

In this work we presented the first extensive statistical analysis of three-dimensional 
simulations of turbulence in collisionless plasma. We studied the case of double- 
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isothermal closure of the CGL-MHD equations in order to compare our results with 
previous isothermal simulations of MHD turbulence. We performed simulations with 
different characteristic sonic and Alfvenic Mach numbers, as well as different pressure 
anisotropics to account for both firehose and mirror instabilities. As main results we 
showed that: 

• we obtained firehose/mirror unstable conditions in all simulations. The unstable 
conditions may be created, even for stable initial conditions, as a consequence of 
the driving and evolutions of turbulent cascade; 

• the supersonic and super- Alfvenic models showed no significant differences when 
compared to standard MffD models. Basically, strong turbulence is able to destroy 
the changes in the topology resulted from the instabilities; 

• the PDF's of density showed broadened profiles for the subsonic and sub-Alfvenic 
cases. The PDF's of velocity showed changes for sub-Alfvenic models. Specifically, 
we obtained an increase on the number high velocity fiows in subsonic models, while 
its decrease for supersonic turbulence; 

• the spectra of density and velocity showed an increase of power in small scales for 
subsonic models; 

• the structure functions of velocity and density fiuctuations revealed that the firehose 
instability tends to isotropize the fluctuations regarding the local reference frame, 
i.e. along the magnetic field lines. On the other hand, the mirror instability 
increases the elongation of the fluctuations along the magnetic fleld lines; 

• the dynamical timescale (r^ ~ SV/L), may also be a good estimate for the 
saturation timescale of the instability growth of most of the wavelengths. This fast 
evolution of the system implies in interesting physical processes in, e.g. interchange 
of energy and acceleration of cosmic rays in the collisionless plasma at intracluster 
medium of galaxies. The growth rate of long wavelengths may be much larger than 
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Appendix A. Linecirization of the Double-Isothermal CGL-MHD Equations 

In this appendix wc derive the dispersion relation for the double-isothermal CGL-MHD 
equations. We start from the isotropic case a| = a\, since it is the case of ideal MHD, 
and then extend the analysis including the anisotropy term by introducing the pressure 
tensor. In this way we can distinguish the role of the anisotropy pressure directly. 
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The ideal double-isothermal CGL-MHD equations with the pressure isotropy 
assumption = a|p = pj_ = a\p results in P = a\pl and can be written in the 
non-conservative form as follows 

% + Vp + pV-t; =0, (A.l) 

at 

p^+ pv-Vv + a\Vp- — {V X B)x B = Q, (A.2) 
at Ait 

dB 

— -Vx{vxB) =0, (A.3) 

where p is the density, v is the velocity field, B is the magnetic field, and a± and ay 
are the sound speeds in the perpendicular and parallel directions with respect to B, 
respectively. 

We assume that all variables can be separated into the uniform and fluctuating 
components, i.e., p — )■ po + 5p, v ^ Vo + 6v, B ^ Bo + 6B. We also assume the lack of 
uniform flow, i.e., vq = 0. Substituting the variable separation in equations (IA.l|) -( IXl3|) 
and removing all non-linear terms leads to the following set of equations 
86 p 



+ po V -dv =0, (A.4) 

Po^ + alV5p - — iVxSB)xBo = 0, (A.5) 

^ - V X (5^ X So) = 0, (A.6) 

Introducing the variable perturbation of the form of 5/(x, t) oc exp [i (k ■ x — ut)] 
results in the above set of equations represented in the Fourier space 

- uj6p + po{k- 6v) =0, (A. 7) 

-ujpo6v+ a]_6pk- — {kx SB) X Bq = 0, (A. 8) 

An 

-ojSB- kx {Svx Bo) =0, (A.9) 

Multiplying equation flA.SP by u, dividing by po and substituting equations flA.7p 
and ( 1A.9I) we obtain the dispersion relation 

- u^6v + ai (fc ■ 6v) k + ^— {[k {Bo ■ Bo) - Bo {Bo ■ k)] {k ■ 6v) 

Anpo 

- k{k-Bo) {Bo ■ 6v) + {k-Bo) {Bo ■ k)6v} = 0, (A. 10) 

Without loosing the generality we can assume that the mean magnetic field is 
parallel to the X direction, i.e. Bq = Box, and k lays in the XY-plane under an angle 6 
with the respect to Bo, i.e. k = k {cos6x + sinOy). We also introduce the Alfven speed 
ca = Bo/y/Anpo- Substituting these assumptions and diving the dispersion relation by 
we can express it in a matrix form 

I -'^ + a\cos^6 a5_ sin 6* cos 6* 
k5v= a\smecose + a\s\i? 6 + c\ | 5vy | = 0. (A.ll) 

V -f^ + cicos2 
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Determinant of the matrix A gives the dispersion relation 



det A = ( -— + c^cos^6' 



- T? l«± + j + Ca«± cos e 



„2 2 



with the eigenvalues 



cos^ 9, 



corresponding to the Alfven wave and 



a]_ + c\± \J (aj_ + c\)^ — Aa']_c\ cos^ 6 



(A.12) 



(A.13) 



(A.14) 



corresponding the fast and slow magnetosonic waves, respectively. 

In the next step we include the pressure anisotropy term in the motion equation 
f lA.2p which for the double-isothermal approximation can be written as 



V 



{p\\ — P±j bb 

„2 „2 



(aj - ai) V ■ [pbb) 



b[b-V)p + ^(b-V)B 



2p 
B 



b-VB] ■ b 



(A.15) 



Substituting the variables separation in this term we obtain 



V ■ 6 {p\i - bb 

(aj - al) Ibo {bo ■ V) (5p + ^ {k .v)5B-'^ [(So ■ V5B) ■ k] , (A.16) 



where bo = Bq/Bo, and the first term corresponds to the perturbation of the density, 
and two remaining terms correspond to the perturbation of magnetic field. Introducing 
the perturbed variables of the form as before, i.e. 5/(x, t) oc exp [i (k • x — cut)], the 
term flAl6|) can be written in Fourier space as follows, 

k ■ S [(j9|| - p^) bb] = (aj - al) {bo ■ k) S^boSp + ^JB - ^ (&o ■ 5B) 6o} • (A.17) 

Next, we multiply it by u and substitute equations (1A.4I) and (1A.6I) to obtain simpler 
relation 



uk ■ 6 



{p\\ - P±) bb = po (aj - a]_) {k ■ 6o)^ 26o (^o • Sv) 



6v 



(A.18) 



Finally, assuming Bq = Bqx and k = k (cos 6x + sin 6y) the term reduces to a very 
simple expression 



--k-5 

P Po 



P\\ 



P±j bb = {26vxbo — Svj (ay — a^) cos^ 6. 



(A.19) 



Rewriting this term in a matrix form we obtain the contribution to the matrix A resulting 
from the presence of pressure anisotropy 



AA 6v 



On — a]_ I COS 








2 „2 



aj ) cos 










(a| — cos' 









SVy 







(A.20) 
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The matrix A + AA takes form 



/ -f;+a}cos2^ 
a\ sin 9 cos 9 

V 



a\ sin 9 cos 9 







Cm COS^ I 









IF + 



— — j cos^ 9 



\ 



, (A.21) 



and its determinant gives the dispersion relation for the CGL-MHD equations 

cos^ d 



detA 



A;2 



+ 



4 



OJ (jJ I 

~ "P" (^^ + «i) + [«! (4 + ai - oy cos^ 9^ - sin^ ^] cos^ ^ |> = 0, (A.22) 



with the eigenvalues 



00 



2 / 2 2 M 2/1 
~ ~ Q± j cos t^, 



(A.23) 



(62 - aj cos2 ^) - a\_ sin^ ^] cos2 9 I , (A. 24) 



corresponding to the Alfven wave and 

corresponding the fast and slow magnetosonic waves, respectively, where h'^ = a\ + c\. 
The derived dispersion relations determine the stability conditions for all characteristic 
waves and growth rates for the firehose and mirror instabilities. 
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